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We provide a fully analytical microscopic theory for the proton correlations in water ice Ih- We 
compute the full diffuse elastic neutron scattering structure factor, which we find to be in excellent 
quantitative agreement with Monte Carlo simulations. It is also in remarkable qualitative agreement 
with experiment, in the absence of any fitting parameters. Our theory thus provides a tractable 
analytical starting point to account for more delicate features of the proton correlations in water 
ice. In addition, it directly determines an effective field theory of water ice as a topological phase. 


The different phases of matter are commonly illus¬ 
trated through the example of ice, water and steam. 
However, common water ice has for a long time been 
known to be a most untypical solid^ exhibiting Paul¬ 
ing’s celebrated ’zero-point entropy^ as the protons in 
fact remaining intricately partially disordered, with only 
the oxygens being associated with a regular arrangement 
on a lattice. 

Nowadays, we can address detailed questions about 
the microscopic nature of water ice. In its most com¬ 
mon form, ice ly, the fourfold coordinated o^gen atoms 
form a hydrogen-bonded wurtzite structure.^ One of the 
Bernal-Fowler ice rules state that each oxygen atom 
has two protons sitting close to it, and the other two 
are further away (Fig. [^. Being subject to these ice 
rules, proton locations are thus not entirely random, 
and the nature of correlations between them has been 
a topic of study for a long time, in theory, simulation 
and experimenfP^H^. In this context, the fact that an 
ice rule can be cast as a conservation law was identified 
as leading to a characteristic pinch-point feature in the 
structure factor of ice.l^lti^ 

However, beyond this long-wavelength insight, little 
analytical progress has occurred, thus providing only lim¬ 
ited analytical backup for the existing extensive numeri¬ 
cal modelling, which does exhibit satisfactory agreement 
with experiment. 

In recent years, much progress in this direction has 
occurred in the context of the family of magnetic com¬ 
pounds known as the spin ice^^^^^, for reviews see 
Refs. I1QI15I These have benefitted much from transfer 
from water ice - including receiving their name. In this 
work, we aim to reverse the direction of this transfer. 
This partially builds on insights gained in considering 
spin ice as a model system for topological phases in con¬ 
densed matter physics - with water ice probably is the 
most ubiquitous material exhibiting topologically non¬ 
trivial behaviour. 

In the following, we first present our analytical theory 
of the proton correlations in water ice, which we con¬ 
trast to simulations and experiment. We then briefly 
put this in the context of topological phenomena in con¬ 
densed matter physics, in particular by deriving an accu- 


FIG. I: Protons in an ice crystal. The locations of the oxygen 
ions approximate a bipartite wurtzite structure (with identi¬ 
cal atoms on the two sublattices, also known as ’hexagonal 
diamond’).^ Like in the case of cubic ice, where they form a 
(cubic) diamond lattice, they are four-fold coordinated with 
the midpoints of the bonds forming tetrahedra. However, 
the tetrahedra are arranged differently: the relation between 
these two lattice types is the same as between hexagonal and 
cubic close packed crystal systems The oxygens are bonded by 
the red hydrogen/deuterium ions. The latter are not ordered 
but are subject to the ice rule that each oxygen has two sitting 
closeby and two further away, so that the {D/H }20 molecules 
‘retain their identity’ even in the solid. This implies an emer¬ 
gent conservation law for the flux denoted by the brown ar¬ 
rows on the midpoints of the bonds, defined to point towards 
the proton on the bond (middle): each oxygen ion sees the 
same amount of flux arriving as leaving. Right: sketch of a 
unit cell with an example of an allowed configuration of the 
eight protons it contains. 


rate value of the coupling constant of the long-wavelength 
action for water ice ly. 

Structure factor: We consider the diffuse scattering 
away from the Bragg peaks, as the peaks are encode the 
‘average’ underlying arrangement of the ions rather than 
the displacement of the protons from it, which is what we 
are interested in. Whereas X-ray scattering is sensitive 
to the charge of the ions and their electrons, and thus is 
not ideally suited for detecting protons, neutrons are a 
suitable probe, especially for deuterated ice, D2O. 

Our theory neglects thermal fluctuations and any 
static disorder that may be present in a given sample and 
focuses on implementing the abovementioned ice rules 
The basic degree of freedom is thus the location Vi^ of 






2 


the proton on bond a of unit cell i. We allow two possi¬ 
ble proton positions for each bond - closer to either one 
or the other oxygen being bonded. This is parametrised 
by an Ising variable S' = ±1: 


where are the midpoints of the oxygen-oxygen bonds 
forming a network of corner-sharing tetrahedra, Si^ are 
the Ising “spins” living at the midpoints of the oxygen- 
oxygen bonds, a is the absolute value of the proton dis¬ 
placement relative to the midpoint, and eo, are unit vec¬ 
tors along the oxygen-oxygen bonds. The index a runs 
from 1 to 8 as the unit cell of the wurtzite structure con¬ 
tains 8 sites. 

The ice rules are enforced by devising a Hamiltonian 
which penalises configurations that do not obey them, 
and studying the resulting ground state correlations. 
How to obtain this has been established in the context 
of spin ic^^. The precise form of the Hamiltonian, with 
details of the analysis sketched below, is relegated to the 
appendix. 

The quantity we compute is the structure factor, which 
is proportional to the neutron scattering cross section: 

iOL 

where q is the wave vector change of the neutron, bia = 
b = const is the scattering length, and the brackets de¬ 
note averaging over all proton configurations obeying the 
ice rules. To express the structure factor in terms of the 



FIG. 2: Theory versus Monte Carlo simulations. Correlation 
function {SnSji) from Monte Carlo simulations (symbols) 
compared to our analytical theory (lines) at zero temperature 
for two inequivalent directions, [100] and [001], and different 
system sizes L. The correlation functions are multiplied by 
the cube of the distance x for better visibility. 

Ising spin correlation function {SiaSji 3 ), we note 

5(q) = 6^ y] 

iaJP 


The spin correlator is absent in one of the following terms 

^giaq(Si„e„-S'j3e^)^ ^ 2(S'jaS'j;3) sin(aqe„) sin(aqe/3) 

+ 2 cos(aqecK) cos(aqe/3) 

whence we finally get the diffuse scattering as 

iS(q) oc ^ (SiaSjp) sm{aqea) 

iajp 

with the computation of the correlation function 
{SiaSj^) carried out in the large-framework (described 
in the appendix). This is based on treating relaxing 
the fixed spin length constraint to one which is only 
obeyed ‘on average’, but which has the advantage of be¬ 
ing exactly soluble. For models closely related to the 
present ones, this has been shown to be an accurate 
appr oximat ion ^ 

Indeed, to demonstrate the accuracy of our evalution 
of the correlation function {SiaSjp), we compare Monte 
Carlo simulations to the analytical expression in Fig. 
We use the notation for a unit cell containing four oxygen 
atoms and eight protons. One finds excellent quantita¬ 
tive agreement in different directions, for different dis¬ 
tances, even capturing finite-size effects faithfully. The 
data for finite-size systems with Nf spins can be obtained 
by explicitly solving for an Nt x dimensional interac¬ 
tion matrix in real space, or equivalently carrying out a 
discrete sum over points in reciprocal space. 

Comparison to experiment: The analytically obtained 
neutron scattering structure factor is compared to neu¬ 
tron resultJ^ in different planes in reciprocal space in 
Figs. The proton displacement was taken to be 
a ~ O.ldSGi^ocP^ where Rqq is the distance between 
the oxygen atoms. 

Analytical results compare reasonably well to the ex¬ 
periment, especially in the first few Brillouin zones. In 
particular, location and orientation of the pinchpoints are 
given correctly, alongside the general structure of regions 
with high and low neutron scattering intensity. 

One important feature of our theory is thus that there 
exists an analytical expression for the correlations which 
can be used as a starting point to understand deviations 
away from the ideal ice model. By providing an ex¬ 
plicit analytical form, this obviates a complex first mod¬ 
elling step (e.g. via numerical fitting procedures to Monte 
Carlo simulations) and therefore allows for a focus on the 
physics beyond the ice rules. 

Such deviations can take many forms. One is a simple 
improvement of our modelling of the scattering form fac¬ 
tor of the proton/deuterium ion, which we have treated 
as a point like scatterer so far, but whose finite extent and 
finite-temperature thermal motion will inevitably lead to 
a change in the structure factor in higher Brillouin zones 
irregardless of any cooperative physics. Similarly, we 
have completely omitted the Bragg peaks due to the oxy¬ 
gen ions, which feature prominently in the experimental 
plots. 
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More interesting are defects in the lattice structure or 
in the bonding - violations of the ice rule - which may be 
static or dynamical. Defects in ice have attracted consid¬ 
erable attention, and much is known about them. In our 
treatment, there is one natural parameter which can be 
added to account for the presence of gauge-charged de¬ 
fects (see below), namely an effective correlation length 
beyond which their presence removes the correlations re¬ 
sponsible for the pinch-points, with the resulting width 
of the pinchpoints reflecting the inverse of this length. 

Long-wavelength theory: The great advantage of our 
method lies in the fact that it works for all distances as 
well as a wide variety of lattices. The peculiar nature of 
the proton correlations become particularly evident from 
the analytics upon closer inspection, with the pinch-point 
features in the structure factor visible even in the absence 
of an explicit identification of an emergent gauge field. 
However, historically, it was noticed already early oil^ 
that the pinchpoints in the neutron scattering were due to 
this conservation law. Basically, the reason this happens 
is because a field obeying only the constraint V • S = 0, 
but which is otherwise free, has no longitudinal degrees of 
freedom, and a pinchpoint is precisely the result of pro¬ 
jecting out the longitudinal component of the (otherwise 
random) proton displacements, in a way which is en tirely 
analogous to how the pinchpoints arise in spin 
[The field B is obtained straightforwardly by drawing an 
arrow from the midpoint of the bond towards the centre 
of charge of the proton. Upon identifying these arrows 
with a unit of an imaginary (lattice) flux, S,the (lattice) 
divergence of this flux vanishes, V • S = 0 - the flux is 
conserved (Fig. [^.] 

Via an analogy to magnetostatics, one obtains an ef¬ 
fective partition function, where the permeability /io of 
the electromagnetic vacuum is replaced by an emergent 
permeability (’stiffness’), K\ 

5 = y . 

This is called topological because the above functional 
goes along with the deconflned ‘Coulomb’ phase of a clas¬ 
sical U(l) gauge theory, which is distinct from a simple 
disordered phase, while not exhibiting any conventional 
(crystalline) order for the protons 

As the wurtzite lattice is not cubic (unlike the case 
of spin ice, which corresponds to ice /c), it is not a pri¬ 
ori obvious that a single stiffness constant is enough to 
describe the action of the gauge field, but it turns out 
that just one stiffness constant is enough, at least ap¬ 
proximately. Indeed, the actual value of this stiffness is a 
priori a free parameter, as it cannot be derived by sym¬ 
metry considerations alone, instead depending on non- 
universal details of the model exhibiting an emergent 
gauge field. Our lattice-based microscopic theory does 
incorporate such information, so that it can be used to 
extract an approximate but accurate estimate from its 
long-wavelength expansion near the pinchpoints. For in¬ 
stance, one can expand Eq. |A1| in the Appendix near 



FIG. 3: Diffuse neutron scattering in water ice. Shown is the 
structure factor in the [hOl] (top), [Okl] (middle) and [hkO] 
(bottom) planes, with dark (light) regions indicating high 
(low) scattering intensity. ^ Lef t: theoretical results for zero 
temperature ^ottom: Eq. |A1| ). Right: experimental result 
from D 2 O ic^ at 20 K (bottom: lOK). Note the pinchpoint 
features at the centres of higher Brillouin zones, where areas 
of high and low scattering intensity meet; these are partially 
obscured by Bragg peaks in the experimental data. 


Qx = Qy = ^ and read off, by direct comparison to the 
predictions from the gauge theory, which have the same 
functional form: 

KRlo « V3/8 . (1) 

This piece of information, together with the general long- 
wavelength form of the theory, is enough to yield the 
full asymptotics of the correlation function in real space. 
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which takes the dipolar form, 


Appendix A: Large-A theory for water ice Ih 




^a^/3 ^(eQjr^ji' ) ) 


47rA 


Note that these results imply that the long-range part 
of the interaction between the protons - which is due to 
the interactions between the dipole moments of the bonds 
arising from the asymmetric location of the protons - 
has an additional component to it which is not of an 
electrostatic origin, but rather due to the contribution 
Eq. This would be present even if the protons were 
electrically neutral, so that neutral particles obeying the 
ice rules would exhibit the same form of the correlations! 

This dichotomy is due to the protons being doubly 
charged-they have an (intrinsic) electric on top of an 
(emergent) gauge charge. The latter associated with the 
emergent conservation law for B. The former is not the 
naive ‘electronic’ charge |e|, but rather related to the 
diverge nce of the electric dipole moment at the location of 
a defectCharged defects in ice are therefore special 
in that they are quasiparticles which not only have an 
irrational electric charge, but also an emergent entropic 
Coulomb charge. 

Final remarks: It is remarkable that a material as well- 
known as water ice should be an embodiment of the topo¬ 
logical physics of unconventional types of order that has 
come into focus relatively recently. For future work, it 
would be most desirable to obtain new neutron scatter¬ 
ing data on a par with the recent X-ray work,l^ in or¬ 
der to allow for a more detailed and quantitative under¬ 
standing of the microscopic proton distribution in water 
ice, beyond the analysis presented there; here we have 
provided a parameter-free analytical theory encoding the 
Bernal-Fowler ice rules, which is applicable at all length- 
scales and can be used as a solid basis for including more 
delicate effects in a simple framework. It is a tantalis¬ 
ing prospect that such analytical theories might be more 
generally useful for modelling diffuse neutron scattering 
experiments. 
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Note added: At the same time as this work, a preprint 
by Benton et al. appeared^ on the proton correlations 
in water ice, in particular in the presence of quantum 
dynamics. 


The (approximate) oxygen positions in water ice can be 
used to define a regular wurtzite structure with a unit cell 
containing four oxygen ions, to which contains 8 protons 
(spins) are associated. The effective interaction matrix 
enforcing the ice rule therefore has dimension 8x8 and 
takes the form 
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where h = cos(g'a,/2\/2), c = cos{qx/^'/2 + 'Jiqy / 4.\/2), 
c = cos{q^/4:V2 - V3qy/4:V2), f = exp[i(ga;/4A/2 + 
gy/4\/6 + qz/2^/3)], g = exp[z(-gj;/4\/2 + qyl^\/^ + 
g^/2\/3)], h = exp[i(-gj^/2\/6 + g2/2\/3)], / = 

exp[i(-ga;/4y2-gj^/4v^+g2/2y3)], g = exp[i(fe/4V2- 
g^/4\/6 + qzl2\pS)\, h = exp[i{qy/2^/6 + g^/2\/5)], * de- 
notes complex conjugation, and T is the temperature. 

The Fourier transform of the Ising spin correlation 
function Gag{x) = (SigSiQ ) can be written in the large- 
N theory approximation 

8 Tj r/t 
n=i ^ ^ 


where Cq^p are the eigenvalues of the interaction matrix 
is a unitary transformation that diagonalizes 
the interaction matrix, and A is determined by the equa¬ 
tion 


8« = E 

q,/i 


1 


where n is the number of unit cells. In the limit of low 
temperature, only the four flat bands of the interaction 
matrix contribute to the correlation function and A —^ 

1 / 2 . 

The structure factor can then be straightforwardly ob¬ 
tained analytically, although it turns out that its form is 
rather lengthy. In some high-symmetry planes, it how¬ 
ever, becomes quite simple. For instance, in the [hkO] 

plane (with Cx = cos and Cy = cos and 

C 5 for cos 'H- sin): 
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cos (\\l\qy) sin - 2^5^ cos sin ( 3 ^!%) + 2CySy sin 


n 2 


3 - cos - 2cos cos (\\[lqy) 
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